Method and Apparatus for Multi-spectral Imaging and Analysis of Skin Lesions and Biological Tissues

ABSTRACT

A multispectral nevoscope that uses specific wavelengths in the visible and infrared spectrum of electromagnetic radiation to transilluminate a skin-lesion or a biological tissue or specimen for imaging and maps multispectral 2-dimensional images into 3-dimensional virtual space for providing 3-D distributions of pre-defined parameters representing the characteristic properties (such as melanin, hemoglobin and deoxyhemoglobin, etc.) of a skin lesion. Methods are disclosed for analyzing and using the characteristic distributions of specific parameters for detection and management of skin-cancers, or characterization of a biological tissue or specimen.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a divisional application of U.S. patent applicationSer. No. 12/539,049 filed Aug. 11, 2009, which claims the benefit ofU.S. Provisional Patent Application No. 61/088,170, filed Aug. 12, 2008,the entireties of which are incorporated herein by reference.

FIELD OF THE INVENTION

This invention relates to methods and apparatus for imaging and analysisof skin lesions and biological tissues.

BACKGROUND OF THE INVENTION

Skin cancer is a significant health problem in the United States. It hasbeen reported that one of five Americans will get some form of skincancer in their lifetime. Currently, nearly half of new cancers arediagnosed as skin cancers. Malignant melanoma, the most fatal skincancer, first forms at the upper layers of skin. When metastasized,cancerous cells from melanoma enter blood vessels and proliferatethroughout the body. Malignant melanoma is highly fatal if not detectedin early stages. However, it can be cured with nearly 100% survival rateif removed at an early stage.

Physicians usually use the “ABCD” rule to determine if a lesion underinvestigation is malignant melanoma. The acronym “ABCD” refers toasymmetry, border, color and diameter, respectively. Malignant melanomatypically has an asymmetrical shape, an uneven border, varied colors anda large diameter. Once a suspicious lesion is excised a diagnosis can beconfirmed by other instruments. However, neither visual inspection usingthe “ABCD” rule nor examination of the excised lesion can provide depthinformation of the skin cancer, which is a crucial signature to gradethe degree of invasion of a skin lesion. Angiogenesis, or increasedblood flow, plays a very important role in detection of melanomas inearly curable stage. Specific patterns of distribution of melanin,oxy-hemoglobin and de-oxy-hemoglobin can lead to characterization ofdysplastic nevi and their potential for transformation into malignantmelanoma in very early phases.

Various light transportation models have been used by researchers toreconstruct information to characterize skin-lesions. For example, aKubelka-Munk model was used to simulate the formation of images ofmelanoma and presented a method to recover blood and melanindistribution in various skin layers. Claridge et al., An inverse methodfor recovery of tissue parameters from colour images, InformationProcessing in Medical Imaging. Springer, Berlin, LNCS2732, pp. 306-317.However, the Kubelka-Munk model is theoretically established in aone-dimensional system with point-based measurements. For more complexgeometries, Monte Carlo simulation or Diffusion Approximation has beenused in optical tomographic modalities for more accuratereconstructions. The commonly adopted strategy for reconstructioninvolves dividing the field of view into a number of voxels and assumingconstant optical properties in each voxel. The optical properties arethen estimated voxel-by-voxel by matching model predicted measurementsto the actual measurements. This is a typically under-determined andill-posed inverse problem as the number of measurements is usually muchless than the number of voxels to be reconstructed. In general, theforward process is a mapping from high dimensional space (unknownoptical properties of voxels) to low dimensional space (limitedmeasurements). Due to the loss of information during the forwardprocess, the solution to the inverse problem is not unique and usuallyhas to be stabilized through various regularization methods. It istherefore difficult to obtain a quantitatively accurate andwell-localized solution. In addition, light photons are quickly diffusedin a turbid medium such as human skin. As a result, there is a strongdependence or similarity between different measurements such thatincreasing the number of measurements would not lead to a dramaticchange in the characteristic behavior of the inverse problem.

In recent years, optical medical modalities have drawn significantattention from researchers. Visible and near-infrared light wavelengthshave been used in surface reflectance, transillumination andtransmission based methods. See, Ganster et al., Computer aidedrecognition of pigmented skin lesions, Melanoma Research, vol. 7 (1997);Seidenari et al., Digital video-microscopy and image analysis withautomatic classification for detection of thin melanomas, MelanomaResearch 9(2), 163-171 (1999); Menzies et al., Automated instrumentationand diagnosis of invasive melanoma, Melanoma Research vol. 7, 13 (1997);Claridge et al., From color to tissue histology: Physics-basedinterpretation of images of pigmented skin lesion, Medical ImageAnalysis, pp. 489-502 (2003); Tomatis et al., Automated melanomadetection: multi-spectral imaging and neural network approach forclassification, Med. Phys. 30(2), pp. 212-221 (2003); Tomatis et al.,Spectro-photo-metric imaging of subcutaneous pigmented lesion:Discriminant analysis, optical properties and histologicalcharacteristics, J. Photochem. Photobiol., B 42, 32-39 (1998). U.S. Pat.No. 5,146,923 discloses a portable nevoscope which provides anoninvasive means to examine a skin lesion in situ, and provides a meansto process and analyze skin lesion data relating to properties such asthickness, color, size, pigmentation, boundary, and texture. Due to thelimited view and limited-angle measurements available via the prior artnevoscope, the intrinsic ill-posed and under-determined nature ofoptical imaging pose problems in reconstructing accurate tomographicinformation.

Consequently there is a need for an improved nevoscope device andmethods of obtaining improved reconstruction results.

SUMMARY OF THE INVENTION

In accordance with various aspects of the present inventionmultispectral imaging systems and methods are provided.

Optical modalities can provide a portable imaging system for routinescreening and monitoring of skin-lesions. Multi-spectral optical imagingusing visible and infrared light wavelengths as disclosed herein canprovide information about physiologically meaningful chromophores suchas melanin, oxyhemoglobin and deoxyhemoglobin through utilization ofdifferences in their wavelength dependent absorption and scatteringcoefficients. The apparatus and methods disclosed herein are generallyapplicable for optical image reconstruction.

In accordance with one embodiment an improved multi-spectral nevoscopeis disclosed providing transillumination for imaging skin lesions fordiagnosing malignant melanoma non-invasively. In one embodiment thedevice comprises substantially a portable optical imaging device thatuses specific wavelengths in the visible and infrared spectrum ofelectromagnetic radiation to transilluminate a skin-lesion or abiological tissue or specimen for imaging and maps multispectral2-dimensional images into 3-dimensional virtual space for providing 3-Ddistributions of pre-defined parameters representing the characteristicproperties (such as melanin, hemoglobin and deoxyhemoglobin, etc.) of askin-lesion. These characteristic distributions of specific parameterscan be analyzed and used for detection and management of skin-cancers,or characterization of a biological tissue or specimen. The deviceallows a background transillumination source for excitation orpreparation of background tissue such as the surrounding skin of a skinlesion or the entire tissue itself for fluoroscopy imaging.

In accordance with an embodiment the device may include multipletransillumination rings for background tissue preparation or excitationfor lesion imaging for optimal penetration and subcutaneous illuminationof skin lesions. The device may include multiple adaptive combinationsof source and receiver channels distributed over the imaging areathrough fiber-optics cables, optical illuminators and filters, andcomputer-controlled image sensors such as CCD arrays. Systems inaccordance with the present invention are of critical value tocharacterize skin lesions and biological tissues for optical and/ormolecular imaging and analyses of associated distributions ofcharacteristic parameters. A series of images obtained with multipleexcitation and source-illumination geometries with multi-spectralfilters may be analyzed by visual inspection/diagnosis and/or 3-Dmapping of distribution of specific parameters such as oxyhemoglobin,deoxyhemoglobin and melanin for diagnostic evaluation andcharacterization of skin lesion or tissue.

The present apparatus may be used for clinical monitoring of skinlesions on patients with high risk of developing malignant melanoma, inaddition to monitoring other skin cancers and conditions including thosedeveloped from allergic reactions in response to drugs, foods and thelike.

In accordance with a further embodiment a shape-based multi-constrainedreconstruction algorithm is disclosed which uses genetic algorithm-basedoptimization methods to find the best possible reconstruction solution.In one embodiment, a skin lesion such as melanoma is modeled as melaninand blood parts, which are delineated by two cubic tensor-productB-spline surfaces. This reduces the number of unknowns to a few controlparameters of the surfaces. The parameters are then coded into a geneticalgorithm to find a solution through global optimization.

In accordance with a further embodiment a multispectral imaging (MSI)method uses plural selected visible wavelengths for transillumination toacquire multiple remittance images. Different wavelengths of light areprojected through fiber-optics-directed ring-light sources fortransilluminating the skin lesion through the surrounding skin area forbackground imaging for calibration, tissue excitation, or tissuepreparation for imaging. The entire remittance signature spectrum usingmultiple light wavelengths improves characterization of skin, dysplasticnevi, melanomas and other skin-lesions. Multiple discrete sources usedin sequential imaging of the skin lesion provide extended image data foruse in characterizing the skin lesion. This characterization may bebased on the visual examination of multispectral transilluminationimages and/or computer based analysis and three-dimensionalreconstruction of the skin lesion.

In accordance with a further embodiment an algorithm is provided inwhich a skin lesion such as melanoma is modeled as melanin, hemoglobinand deoxyhemoglobin.

BRIEF DESCRIPTION OF THE DRAWINGS

To assist those of ordinary skill in the relevant art in making andusing the subject matter hereof, reference is made to the appendeddrawings, wherein:

FIG. 1 is a schematic diagram of a prior art nevoscope apparatus;

FIG. 2 is a schematic diagram of an illumination and imaging system inaccordance with at least one embodiment of the present invention;

FIG. 2A is a schematic diagram of a layout of contiguous fiber-opticsbundles-based imaging geometry, including fiber bundles with receiverchannels and illuminating channels with N-to-1 geometry, in accordancewith at least one embodiment of the present invention;

FIG. 3 is a schematic diagram of a sensor face plate in accordance withat least one embodiment of the present invention;

FIG. 4A is a schematic of an epi-illumination mode for the apparatus ofFIG. 1;

FIG. 4B is a schematic of a trans-illumination imaging mode for theapparatus of FIG. 1;

FIG. 5 is a schematic diagram of a method of shape-based reconstructionin accordance with at least one embodiment of the present invention;

FIG. 6A is a schematic diagram of a discretization strategy of detectorspace in accordance with at least one embodiment of the presentinvention;

FIG. 6B is a schematic diagram of a discretization strategy of theinterrogated tissue medium in accordance with at least one embodiment ofthe present invention;

FIG. 7 is a graphical representation of a shape-based model of malignantmelanoma in accordance with at least one embodiment of the presentinvention;

FIGS. 8A-8E are graphical depictions of reconstruction results: FIG. 8Areflects double-surface model results (left: first surface, right:second surface) and FIGS. 8B-8E reflect reconstructed surfaces withdifferent constraints in accordance with at least one embodiment of thepresent invention;

FIG. 9 is a graphical depiction of convergence analysis of the geneticalgorithm with various constraints in accordance with at least oneembodiment of the present invention;

FIGS. 10A and 10B are graphical depictions of a deformation processduring optimization in accordance with at least one embodiment of thepresent invention (from left to right rowwise and top to bottomcolumn-wise): FIG. 10A depicts the first surface and FIG. 10B depictsthe second surface;

FIG. 11 is a schematic diagram of a 3-dimensional feature reconstructionfrom multispectral images in accordance with at least one embodiment ofthe present invention;

FIG. 12 is a schematic diagram of a computer classification system usingadaptive fuzzy clustering and portioning in accordance with at least oneembodiment of the present invention;

FIG. 13A is a schematic diagram of partitioning of the feature spacewith hyperplanes in accordance with at least one embodiment of thepresent invention; and

FIG. 13B is a schematic diagram of a winner-takes-all strategy usingfuzzy partitions in accordance with at least one embodiment of thepresent invention.

It should be noted that the appended drawings illustrate only typicalembodiments of this invention and are therefore not to be construed aslimiting of its scope, for the invention may admit to other equallyeffective embodiments. Where possible, identical reference numerals havebeen inserted in the figures to denote identical elements.

DETAILED DESCRIPTION OF THE INVENTION

In the following description, for purposes of explanation, specificnumbers, materials and configurations are set forth in order to providea thorough understanding of the invention. It will be apparent, however,to one having ordinary skill in the art that the invention may bepracticed without these specific details. In some instances, well-knownfeatures may be omitted or simplified so as not to obscure the presentinvention. Furthermore, reference in the specification to phrases suchas “one embodiment” or “an embodiment” means that a particular feature,structure or characteristic described in connection with the embodimentis included in at least one embodiment of the invention. The appearancesof phrases such as “in one embodiment” in various places in thespecification are not necessarily all referring to the same embodiment.

Now referring to FIG. 1 a schematic of a prior art Nevoscope apparatus 2employing white light-based transillumination is shown including anoptical lens 10, lens mount bracket 20, surface light fibers 30,transillumination fibers 40, a ring light 50 and monocoil 60. As shownthe nevoscope is positioned over skin 100.

Now referring to FIG. 2, an embodiment of an improved Nevoscopeapparatus 200 is shown. Apparatus 200 includes a housing 210, face plate220, transillumination outer ring fiber cable 240, main imaging areafiber cables 250 and 252, focusing lens/polarizer 260, image sensor CCD270, source mask 280 and multispectral filter 290. Source mask 280 andmultispectral filter 290 may be rotatable and/or selectable, such as bymicroprocessor-based computer control. Apparatus 200 may include atransillumination mask and/or filter which may be microprocessor-basedcomputer controlled. A face plate height adjuster 230 such as forepi-illumination may be included in the apparatus 200.

Now referring to FIG. 2A, contiguous fiber-optic bundles with receiverchannels 306 and illuminating channels 308 with N-to-1 geometry areemployed in an alternative nevoscope apparatus. Instead of anincandescent illuminator with a set of filters to produce visible andnear-infrared wavelengths, multiple multispectral surface mount lightemitting diodes (LED) (such as are commercially available from LumexInc.) with an appropriate multiplexed LED driver (such as LTC3219 fromLinear Technology) to drive the LEDs can be used. A square pulse signalmay be used to turn on the LEDs. A control signal can be used to directthe turn-on time (preferably in the order of a few milliseconds) to aparticular LED at a time in a sequential manner. The control and theturn-on signal pulses can be synchronized to the camera frame rate inorder to ensure each picture frame corresponds to a single LEDillumination.

Fiber optic bundles (such as those made from high quality silica withbetter than 99% transmission from 450 nm to 960 nm, available throughSCHOTT North America) can be used also in a contiguous manner ratherthan distributed manner for illumination and receiver channels in aparticular geometry (such as in alternate mode or N-to-1 mode in whichthere are N number of receiver channels for each illumination channel;where N is positive integer preferably N=1,2, . . . 64; when N=1, itbecomes the alternate mode).

Fiber optic bundles can be directly divided into illumination andreceiver multi-core channels where multi-core illumination channels areconnected with multispectral LEDs through multiplexed LED drivers(preferably LTC3219, Linear Technology) and receiver channels areconnected to a CCD camera (such as Sony ICX415 CCD) through a focusingoptical lens.

Now referring to FIG. 3, face plate 220 includes outer ring illuminators222 and 224 and imaging area 226 with distributed sources in alternatepositions in a matrix corresponding to the source mask. Outer ringilluminator 222 includes shorter wavelength transillumination fiberchannels 223 oriented at 45 degree convergent beam for multispectral orspecific wavelength-based background image excitation. Outer ringilluminator 224 includes longer (relative to illuminator 222) wavelengthtransillumination fiber channels 225 oriented at 45 degree convergentbeam for multispectral or specific wavelength-based background imageexcitation. Imaging area 226 includes plural illumination fiber channels227 operably connected to an illuminator-filter assembly and pluralreceiver fiber channels 228 operably connected to the CCD image sensor270.

Transillumination imaging is achieved through 45 degree convergent beamsthrough fibers distributed along outer ring(s) 222 and 224 as describedabove with separate illumination, transilumination masks 280 andmultispectral filter 290 selection and control. The describedtransillumination using any selected optical wavelength can be used fortissue excitation such as for fluoroscopy and/or simpletransillumination imaging of background skin/medium.

Face plate 220 may be removably attached to housing 210 and is operablyconnected to a removable lens/polarizer 260 that provides an interfaceto the imaging areas. The fibers 227 and 228 are distributed over theimaging area in a matrix that can be radially symmetric or rectangular(as shown). The fibers 227 and 228 from the imaging area 226 split intofiber cables 250 and 252 with the same access to imaging area 226.

Now referring to FIG. 2, for imaging a specific imaging mask may becreated and used in the illumination path of fiber cable 252. The maskhas openings for the desired fibers to be used as source locations tosend light into the specimen/skin-tissue 100. As will be apparent to theskilled artisan, many possible schemes for illumination may be employedby fiber placement.

For multi-spectral imaging, any optical filter of a specific wavelengthpass or band filter can be used in the illumination pathway of fibercable 252 and selected through a computer controlled interface.

For recording images, an appropriate mask is used to receive the lightfrom the fibers that are not used for illumination. The received lightis passed through a focusing lens/polarizer 260 to a CCD sensor 270 toform images and record measurements. The corresponding receiver mask canbe appropriately selected through a computer controlled interface.

Now referring to FIG. 4A, in an epi-illumination mode, light beams L arereflected from above the skin surface 100 and diffused reflected lightRL is collected by CCD sensor 270 (not shown) through the opticalassembly of cross-polarizer and magnifying lens. This image carriesapparent characteristic of a lesion that can be used for automaticdiagnostic algorithm in terms of the “ABCD” rule. Now referring to FIG.4B, in a trans-illumination mode, photons of white-light spectrum or aspecific wavelength are directed by a transilluminator ring light suchas described in FIG. 3 providing an optical interface for light L toenter into the surrounding area of a skin-lesion with a cone-beam makinga forty-five degree angle with respect to the normal of skin 100. Theback-scattered diffused light DL re-emerges from the skin and capturedby the CCD sensor 270 (not shown) such as a camera through the opticalassembly. This image contains the information about absorption andscattering properties of the chromophores of underlying skin layers andlesion. Within the optical tomography framework, it is possible toretrieve the distribution of melanin and blood as two key signaturevariable for early detection of malignant melanoma. In one embodiment, ashape-based multi-constrained algorithm is applied for reconstructionresults, which algorithm overcomes the intrinsic problems associatedwith reconstructing accurate tomographic information.

Now referring to FIG. 5 a method of reconstruction is described in termsof Nevoscope transillumination images. The method minimizes thedifference between the actual (ΔM^(real)) and predicted (ΔW^(cal));where ΔM^(real) is the actual measurement vector obtained frommultispectral images and ΔM^(cal) is the computed measurement vectorobtained from the reconstructed images. measurement by employing geneticalgorithms providing global searching. A linearized forward model isadopted and evaluated by Monte Carlo simulation in terms of typicaloptical properties of normal skin. In one aspect, malignant melanoma isrepresented by shapes of its melanin part and blood part. Theseparameters are grouped into genetic algorithms.

To develop a reconstruction strategy, a forward model is required torelate the measurement to the optical properties of tissue underinvestigation. Regardless of what kind of imaging geometry is used, anoptical system may be described as

M=F(x)  (1)

where M is the measurement and F is a forward model. x is a distributionof unknown optical properties. Given a reasonable initial guess x₀ ofthe background optical properties, Equation (1) may be expanded into

M=F(x ₀)F′+(x x ₀)½+F″(x ₀)(x x ₀)  (2)

where F′ and F″ are first order and second order Frechet derivativesrespectively.

Let ΔM=M−F(x₀) and Δx=x-x₀ , Equation (2) may be re-arranged as

ΔM=F′Δx½F″Δx  (3)

The discrete form of Equation (3) turns out to be

Δ{right arrow over (M)}=J Δ{right arrow over (x)}+½H Δ{right arrow over(x)}+  (4)

Here, J is the Jacobian matrix and H is the Hessian Matrix. Δ{rightarrow over (M)} is the measurement vector and Δ{right arrow over (k)} isthe vector that gives the variations from the background {right arrowover (x)}₀ . Neglecting higher order terms in Equation (4), a simplifiedlinear system can be expressed as

Δ{right arrow over (M)}=Δ{right arrow over (x)}  (5)

The formulation in terms of Equation (5) leads to linear opticaltomography which is also known as “difference imaging” with twomeasurements. One is for background tissue (that is, x₀) and another isfor abnormal tissue (that is, unknown x). The difference is then fed tothe reconstruction algorithm to obtain the optical properties. In thismethod, such a linear approach is adopted for nevoscope and the Jacobianmatrix is extracted by Monte Carlo simulation.

Monte Carlo (MC) simulation is viewed as the “gold standard” to modelinglight transportation and has been extensively used by researchers inoptical communities. One of the major types of MC simulations predictsthe trajectories of photons in an absorption-free medium then attenuatesthe weights of photons in terms of microscopic Beer's law. Using thisapproach, the trajectories of photons can be stored and used repeatedlywhich speeds up the calculation and absorption properties in the modeland can be easily updated. A Monte Carlo procedure based on thisapproach and tuned to Nevoscope geometry is described as follows.

A nevoscope uses a ring light with radius r and incident angle 45 degreewith respect to the normal of skin surface. In the coordinate systemused in the experiments that follow, skin surface is represented by x-yplane at z=0 and z axis points inward into the medium. The light sourceis realized by launching a photon at a random position (x₀, y₀, 0) onthe ring, where

x ₀ =r cos(2πξ) and y _(o) =r sin(2πξ)  (6)

where, ξ is a uniform random number from 0 to 1. Because the discrepancyof refractive indices between air and skin, the photon will change itsincident direction which is governed by Snell's law

n _(i) sin(α_(i))=n _(t)sin(α_(t))  (7)

As the photon loses part of its energy due to specular reflection,considering each photon has normalized energy 1, the remaining energy isgiven by

$\begin{matrix}{w_{0} = {1 - {\frac{1}{2}\left\lbrack {\frac{\sin^{2}\left( {\alpha_{i} - \alpha_{t}} \right)}{\sin^{2}\left( {\alpha_{i} + \alpha_{t}} \right)} + \frac{\tan^{2}\left( {\alpha_{i} - \alpha_{t}} \right)}{\tan^{2}\left( {\alpha_{i} + \alpha_{t}} \right)}} \right\rbrack}}} & (8)\end{matrix}$

For Equations (7) and (8), a_(t)=45° is the incident angle and α_(t) isthe transmitted angle. =1 is the refractive index of the ambient and n,is that of stratum corneum.

After the photon is launched, it undergoes multiple scattering beforetermination. The probability of path length that a photon wouldencounter a scattering event is given by Sobol, The Monte Carlo method(Chicago, Ill.: University of Chicago Press), 1974:

p(l)=μ_(s)exp(μ_(s) l)−exp(μ_(α) l)−μ_(α)+exp(μ_(α) l)−exp(μ_(s) l)  (9)

To predict the photon trajectory first, absorption may be neglectedtemporarily. That is, the simulation is done with an absorption freeskin model. Thus Equation (9) reduces to

p(l)=μ, exp(μ_(s) l)  (10)

So the path length between scattering events is derived as

$\begin{matrix}{l = \frac{\ln (\xi)}{\mu_{s}}} & (11)\end{matrix}$

When a photon completes its one step, it will scatter into a newdirection. The direction is characterized by deflection and azimuthalangles. The azimuthal angle is sampled from a uniform distribution from0 to 2π but the deflection angle has to be sampled from a phase functionwhich represents the typical forward scattering of biological tissue. Aphase function may be expressed as

$\begin{matrix}{{p\left( {\cos \mspace{11mu} \theta} \right)} = {\frac{1}{4\pi}\frac{1 - g^{2}}{\left( {1 + g^{2} - {2g\mspace{11mu} \cos \mspace{11mu} \theta}} \right)^{\frac{3}{2}}}}} & (12)\end{matrix}$

In addition, it is quite likely that a photon will travel acrossinternal refractive index mismatched interfaces or hit the skin surfaceduring its transportation. If a photon hits the surface, its energy canbe separated into two parts: one part escapes the surface and the otherpart is internally reflected to propagate with reduced energy. The ratiobetween internal reflected energy and the original energy is given byFresnel reflection coefficient.

$\begin{matrix}{{R\left( \alpha_{i} \right)} = \left\{ \begin{matrix}\frac{\left( {n_{t} - n_{i}} \right)^{2}}{\left( {n_{t} + n_{i}} \right)^{2}} & \left( {{{if}\mspace{14mu} \alpha_{i}} = 0} \right) \\{\frac{1}{2}\left\lbrack {\frac{\sin^{2}\left( {\alpha_{i} - \alpha_{t}} \right)}{\sin^{2}\left( {\alpha_{i} + \alpha_{t}} \right)} + \frac{\tan^{2}\left( {\alpha_{i} - \alpha_{t}} \right)}{\tan^{2}\left( {\alpha_{i} + \alpha_{t}} \right)}} \right\rbrack} & \left( {{{if}\mspace{14mu} 0} < a_{i} < {\sin^{- 1}\left( \frac{n_{i}}{n_{t}} \right)}} \right) \\1 & \left( {{{if}\mspace{14mu} {\sin^{- 1}\left( \frac{n_{i}}{n_{t}} \right)}} < a_{i} < \frac{\pi}{2}} \right)\end{matrix} \right.} & (13)\end{matrix}$

In Equation (13), a, and a, are angle of incidence and transmittancerelated by Snell's law. n_(i) and n_(t) are medium of incidence andtransmittance respectively.

In the case of a photon crossing internal interfaces, the “Russianroulette” approach proposed by Wang et al., MCML-Monte Carlo modeling oflight transport in multi-layered tissues, Computer Methods and Programsin Biomedicine 47, 131-146 (1995), may be employed. Once the photon hitsan interface, the Fresnel reflection coefficient R(α_(i)) is calculatedand compared to a uniform-distributed random number ξ. If R(α_(i))>ξ,the photon will cross the boundary. Otherwise, it is totally reflected.

Following the rules described above, a photon is launched from the ringsource and propagates through the optical skin model until terminated ordetected. If its total path length exceeds some preset value, the photongets terminated. If it is captured by a detector, its energy isre-weighted in terms of its trajectory as

$\begin{matrix}{{w_{d} = {\left( {1\mspace{11mu} - {R\left( \alpha_{h} \right)}} \right)w_{0}{\sum\limits_{i = 1}^{h - 1}\; {R\left( \alpha_{i} \right)}}}},} & (14)\end{matrix}$

where it is assumed the photon hits the skin surface h times during itstransportation.

After all trajectories of photons are recorded for a specific detector,the total received intensity is therefore given in terms of microscopicBeer's law

$\begin{matrix}{M = {\sum\limits_{i = 1}^{p}\; {w_{d}^{i}{\sum\limits_{j = 1}^{q_{p}}\; ^{{- {\mu_{a}{(j)}}}l_{j}^{i}}}}}} & (15)\end{matrix}$

where w_(d) ^(i) of i th photon is given by equation (14), μ_(a) (j) isthe absorption coefficient along the path l/of ith photon. It is assumedthat the detector received p photons while each of them used q_(p)steps. The wavelength dependent absorption coefficients are obtainedfrom published data. See, R. Anderson et al., The optics of human skin,Journal of Investigative Dermatology, Vol. 77(1), 13-19 (1981); A.Bashkatov et al., Optical properties of human skin, subcutaneous andmucous tissues in the wavelength range from 400 to 2000 nm, Journal ofPhysics D: Applied Physics, Vol. 38, 2543-2555 (2005); M. Van Gemert etal., Skin Optics, IEEE Transactions on Biomedical Engineering, Vol.36(12), 1146-1155 (1989). For purposes of the present experiments aseven-layered skin model is employed which includes stratum corneum (20μm), epidermis (80 μm), dermis (150 μm), upper blood net dermis (80 μm),reticular dermis (1500 μm), deep blood net dermis (100 μm) andsubcutaneous fat (6000 μm). See, I. Meglinski et al., Quantitativeassessment of skin layers absorption and skin reflectance spectrasimulation in the visible and near-infrared spectral regions, Physiol.Meas, 23, 741-753 (2002).

The absorption coefficient of each layer is a combined effectcontributed by several absorbers such as blood, melanin, water and thebaseline skin which is free of these absorbers. The absorptioncoefficient of baseline skin is given by

μ_(α) ^(baseline)(λ)=7.84×10⁸ ×l ^(−3.255) (cm⁻¹)  (A.1)

Here, λ is the wavelength measured in nanometers. The absorption ofblood is determined by its oxygenation and fraction of red cells inblood. For example, we assume that haematocrit is Ht % and S %hemoglobin is oxy-hemoglobin. Then the absorption of blood is given as

μ_(α) ^(blood)(λ)=Ht×(1-S)×μ₆₀ ^(Hb)(λ)+Ht×S×μ ₆₀ ^(HbO) ₂(λ)  (A.2)

Here μ₆₀ ^(Hb)(λ) and μ_(α) ^(HbO) ₂(80 ) are absorption coefficients ofoxy-hemoglobin and deoxy-hemoglobin respectively. The absorption spectracome from compiled data. See, Prahl, Oregon Medical Laser Center,Optical absorption of hemoglobin,http://omlc.ogi.edu/spectra/hemoglobin/index.html (2006).

The haematocrit has a typical value 45% . The absorption of melanin isgiven by Meglinski et al.

μ_(α) ^(melanin)(80 )=5×10¹⁰λ^(3.33)(cm⁻¹)  (A.3)

In healthy skin, melanin only exists in epidermis layer. Its fractionranges from 1% to 3% for light skinned Caucasians, from 11% to 16% forwell-tanned Caucasians and Mediterraneans, and from 18% to 43% fordarkly pigmented Africans. The presence of melanin in deeper skin layerslike dermis usually indicates the malignancy of melanoma. The absorptionof water is known. Absorption coefficients of two blood free layersstratum corneum and epidermis are predicted by

μ_(α) ^(stratum)(λ)=((0.1-0.3×10⁻⁴λ)+0.125μ_(α) ^(baseline)(λ))(1-C_(H)₂ _(O)) +C_(H) ₂ _(Oμ) _(α) ^(H) ² ^(O)(λ)  (A.4)

μ_(α) ^(epidermis)(λ)=(C_(melanin)μ_(α) ^(melanin)(λ)+C_(H) ₂ _(O)μ_(α)^(H) ² ^(O)(λ) +(1-C_(melanin)-C_(H) ₂ _(O)) μ_(α)^(baseline)(λ))  (A.5)

In the above equations, (C_(melanin) is the fraction of melanin inepidermis and C_(H) ₂ _(O) is the fraction of water in a specific layer.Absorption coefficients of four dermis layers and subcutaneous fat aredescribed by one formula (A.6) while their values depend on how muchblood and water they have.

μ_(α) ^(bloodlayer)(λ)=C_(blood)μ_(α) ^(blood)(λ)+C_(H) ₂ _(O)μ_(α) ^(H)² ^(O)+(1-C_(blood)-C_(H) ₂ _(O))μ₆₀ ^(baseline)  (A.6)

Additional optical properties of skin include scattering coefficients,anisotropy factor and refractive index. Representative values of theseoptical properties and fractions of absorbers in each skin layer arereported by Meglinski, supra.

Now referring to FIG. 6A, in accordance with the experiments, thedetector space of the nevoscope is divided into k number of rings 202(for example, three rings 202 as shown in the FIG. 6A) with equal width.The innermost circle 204 contains four elements or detectors 206 withequal size and other rings 202 are split into a number of elements whichmaintain the same area as the ones in the innermost circle 204. Usingthis discretization strategy makes use of the symmetrical geometry of anevoscope. When evaluating Jacobian matrix, instead of obtaining one MCsimulation for each detector, since MC is extremely time-consuming, andnoting that all detectors 206 on the same ring actually have similartrajectories, the trajectories for one detector 206 on that ring 202 isrecorded and trajectories of other detectors 206 are generated byrotating recorded trajectories. So eventually, only k independentsimulations are required to evaluate Jacobian matrix.

Resolution of the reconstructed image is dependent on the size of theJacobian matrix. Now referring to FIG. 6B, in the experiments, theinterrogated volume is divided into N×N×M voxels 300. Each voxel 400 hasthe size Δx×Δy×Δz and has homogeneous absorption coefficient. Where thetotal number of detectors is T, the measurements of background normalskin is then given by

$\begin{matrix}{{{{M(s)} = {{\sum\limits_{i - 1}^{p_{s}}{w_{skin}^{i}{\sum\limits_{i = 1}^{p_{s}}{\left( {w_{d}^{i} = {\sum\limits_{j = 1}^{N \times N \times M}\; ^{{- {\mu_{a}{(j)}}}l_{j}^{i}}}} \right)\mspace{14mu} s}}}} = 1}},{\ldots \; T}}\mspace{11mu}} & (16)\end{matrix}$

where p_(s) is total number of photons received by the detector s. l_(j)^(i) is the path length of photon i in the voxel j. μ_(α)(j) is theabsorption coefficient of voxel j and is extracted from opticalproperties of normal skin.

For difference imaging, the Jacobian matrix is given by

$\begin{matrix}{J = {\begin{bmatrix}\frac{{\partial\Delta}\; {M(1)}}{\partial{{\Delta\mu}_{a}(1)}} & \frac{{\partial\Delta}\; {M(1)}}{\partial{{\Delta\mu}_{a}(2)}} & \; & \; & \frac{{\partial\Delta}\; {M(1)}}{\partial{{\Delta\mu}_{a}\left( {N \times N \times M} \right)}} \\\bullet & \bullet & \; & \; & \bullet \\{\bullet \;} & \; & \bullet & \; & \bullet \\\bullet & \; & \; & \bullet & \bullet \\\frac{{\partial\Delta}\; {M(T)}}{\partial{{\Delta\mu}_{a}(1)}} & \frac{{\partial\Delta}\; {M(T)}}{\partial{{\Delta\mu}_{a}(2)}} & \; & \; & \frac{{\partial\Delta}\; {M(T)}}{\partial{{\Delta\mu}_{a}\left( {N \times N \times M} \right)}}\end{bmatrix}_{{{\Delta\mu}_{a}{(j)}} = 0}}} & (17)\end{matrix}$

Each term of Equation (17) is evaluated by recorded photon history. Whena tumor grows on the normal skin, variation in measurements ΔM(s) may beexpressed as

$\begin{matrix}{{{\Delta \; {M(s)}} = {\sum\limits_{i = 1}^{p_{s}}\left( {w_{skin}^{i}{\sum\limits_{j = 1}^{N \times N \times M}\; ^{{- {{\Delta\mu}_{a}{(j)}}}l_{j}^{i}}}} \right)}}\mspace{25mu}} & (18)\end{matrix}$

Taking the derivative with respect to Δμ_(α)(j) and set Δμ_(α)(j) tozero results in

$\begin{matrix}{{\frac{{\partial\Delta}\; {M(s)}}{\partial{{\Delta\mu}_{a}(j)}}_{{{\Delta\mu}_{a}(j)} = 0}} = {\sum\limits_{i = 1}^{p_{s}}\left( {w_{skin}^{i}l_{j}^{i}} \right)}} & (19)\end{matrix}$

The Jacobian matrix used in the experiments is evaluated on the normalskin model introduced hereinbelow with k=20, N=32 and M=20 which dividesthe nevoscope detector plane of 1.2 cm diameter into 1588 elements anddivides the 1.2 cm×1.2 cm ×2000 μm field of view into 32×32×20 voxels.

Inverse Problem

The Equation (5) can be solved by Singular Value Decomposition (SVD),Algebraic Reconstruction Technique (ART), Conjugate Gradient method (CG)and their variants. However, due to the ill-posed and under-determinednature of the inverse problem, it is unlikely to reconstructquantitatively accurate and well-localized absorption coefficients.Therefore, shape-based optical tomography is preferred, wherein theabnormalities are assumed to have piece-wise constant optical propertiesin the compact supported regions and the boundaries of these regions canbe effectively approximated by some shape functions with the limitednumbers of control parameters. As to the reconstruction problem, theunknowns just include a few numbers of control parameters and,sometimes, the piece-wise constant optical properties.

Shape Representation of Malignant Melanoma

A melanoma is essentially a three-dimensional object with distributedpigmentation and chromophores. Malignant melanoma is a result ofuncontrolled replication of melasome cells sitting on the basal layer ofepidermis. The shape of the melanoma is hence bounded by the epidermislayer. What is needed to describe the lesion is therefore reduced to 2Dsurfaces in the 3D domain.

Breaking the malignant melanoma into parts, such as the melanin and theblood parts, permits representation of the malignant melanoma with 2Dsurfaces. The melanin part is a 3D region bounded by a single surfaceand the epidermis layer. Within the region, the optical properties areconstant and the only absorber is melanin. A second surface which sitsbelow the first surface is used to represent the blood part. The regionbounded by the first surface and the second surface is blood only. Thismodel mimics a skin-lesion and its x-z intersection is shown in FIG. 7.

The first surface is represented as ƒ₁(x, y) which corresponds to thedepth of lesion from the epidermis layer at the position (x, y). Thecontinuous surface is represented with limited parameters. First, a N×Nrectangular grid is positioned to lie over the epidermal layer. Second,the function ƒ₁(x, y) is sampled to N×N discrete values ƒ_(d1)(X, Y).Here, (x, y) is continuous and (X, Y) is N×N numbers of discretesampling positions. Third, the discrete values are interpolated by thecubic tensor-product B-spline which satisfies the following condition:

$\begin{matrix}{{f_{d\; 1}\left( {X,Y} \right)} = {\sum\limits_{i = 1}^{N}\; {\sum\limits_{j = 1}^{N}\; {{c_{1}\left( {,j} \right)}{\beta^{3}\left( {{X - },{Y - j}} \right)}}}}} & (20)\end{matrix}$

and the original function ƒ₁(x, y) can then be approximated

$\begin{matrix}{{{by}\mspace{14mu} {f_{B\; 1}\left( {x,y} \right)}} = {\sum\limits_{i = 1}^{N}\; {\sum\limits_{j = 1}^{N}\; {{c_{1}\left( {,j} \right)}{\beta^{3}\left( {{x - },{y - j}} \right)}}}}} & (21)\end{matrix}$

where

β³(x−i, y−j)=β³(x−i)□β³(y−j)  (22)

is the tensor product of one-dimensional cubic B-spline basis β³(x−i)and β³(y−j), and c₁(i, j) is B-spline coefficient.

The basis function of B-spline comes from repeated convolution of thebox function. It has compact support so that the computation burden isalleviated. The cubic B-spline is continuous up to second derivativewhich leads to the smoothness of the interpolated surface and in essenceputs regularization effect on the reconstructed surface in ourshape-based approach. In the end, for the melanin part, the parametersto determine its shape are reduced to N×N discrete values ƒ_(d1)(X, Y)or, equivalently, the corresponding B-spline coefficients c₁(i, j).

Similarly, the second surface can be defined by N×N discrete valuesƒ_(d2)(X, Y) or, equivalently, z_(d l (X ,Y)=ƒ) _(d2)(X,Y) ƒ_(d1)(X,Y)-which is the thickness of blood region between the first surface andthe second surface.

Shape-Based Reconstruction

To reconstruct the surfaces and piecewise constant optical properties,the continuous surface representation should be incorporated into theforward photon transportation model. In our study, the linearizedforward model is kept intact but the continuous representation issampled into the discrete vector of unknowns Δ{right arrow over (r)}.That is,

Δ{right arrow over (M)}=J□Δ{right arrow over (x)}=J□S(ƒ_(d1)(X, Y), mp1,z _(d)(X,Y),bp2)  (23)

where S(□) is the sampling function which converts the continuous shaperepresentation to the voxel-based optical properties in forward model,mp1 is the fraction of melanin and bp2 is fraction of blood. The inverseproblem can therefore be formulated as minimizing the objective function

F _(obj)=1 ₂∥Δ{right arrow over (M)}−J□Δ{right arrow over (x)}∥²=1₂∥Δ{right arrow over (m)}−J□S(ƒ_(d1)(X,Y), mp1,z _(d)(X,Y),bp2)∥²  (24)

The unknowns of this inverse problem are reduced to ƒ_(d1)(x, y),z_(d)(X,Y), mp1 and bp2. The multi-spectral shape reconstruction using Nwavelengths can be formulated as a multi-objective optimization problemand its objective function is given as

F _(obj)=α₁ □F _(obj) ^(λ) ¹ α₂ −□F _(obj) ^(λ) ² ... +α_(N) −□F _(obj)^(λ) ^(N)   (25)

where, {60 ₁, α₂, . . . α_(N)} is a set of coefficients to balance thecontributions from different single-wavelength objective functions.

The wavelength dependent objective functions are given as

$\begin{matrix}{{F_{obj}^{\lambda_{1}} = {{\frac{1}{2}{{{\Delta \; {\overset{\rightharpoonup}{M}}^{\lambda_{1}}} - {J^{\lambda_{1}}\Delta {\overset{\rightharpoonup}{\; x}}^{\lambda_{1}}}}}^{2}} = {\frac{1}{2}{{{\Delta \; {\overset{\rightharpoonup}{M}}^{\lambda_{1}}} - {J^{\lambda_{1}}\bullet \; {S^{\lambda_{1}}\left( {{f_{d\; 1}\left( {X,Y} \right)},{{mp}\; 1},{z_{d}\left( {X,Y} \right)},{{bp}\; 2}} \right)}}}}^{2}}}}{F_{obj}^{\lambda_{2}} = {{\frac{1}{2}{{{\Delta \; {\overset{\rightharpoonup}{M}}^{\lambda_{2}}} - {J^{\lambda_{2}}\Delta {\overset{\rightharpoonup}{\; x}}^{\lambda_{2}}}}}^{2}} = {\frac{1}{2}{{{\Delta \; {\overset{\rightharpoonup}{M}}^{\lambda_{2}}} - {J^{\lambda_{2}}\bullet \; {S^{\lambda_{2}}\left( {{f_{d\; 1}\left( {X,Y} \right)},{{mp}\; 1},{z_{d}\left( {X,Y} \right)},{{bp}\; 2}} \right)}}}}^{2}}}}\mspace{20mu} \bullet \mspace{20mu} \bullet \mspace{20mu} \bullet {F_{obj}^{\lambda_{N}} = {{\frac{1}{2}{{{\Delta \; {\overset{\rightharpoonup}{M}}^{\lambda_{N}}} - {J^{\lambda_{N}}\Delta {\overset{\rightharpoonup}{\; x}}^{\lambda_{N}}}}}^{2}} = {\frac{1}{2}{{{\Delta \; {\overset{\rightharpoonup}{M}}^{\lambda_{N}}} - {J^{\lambda_{N}}\bullet \; {S^{\lambda_{N}}\left( {{f_{d\; 1}\left( {X,Y} \right)},{{mp}\; 1},{z_{d}\left( {X,Y} \right)},{{bp}\; 2}} \right)}}}}^{2}}}}} & (26)\end{matrix}$

Various optimization techniques have been reported to solve theshape-based reconstruction problems but all have drawbacks. Hence, agenetic algorithm is used because the gradient need not be evaluatedwhich simplifies the computation and provides reliability. Geneticalgorithm is one of the most popular methods used in seeking globalminimal. Among the global optimization techniques, genetic algorithmprovides a reasonable convergence rate due to its implicit parallelcomputation. See, Whitely, A Genetic Algorithm Tutorial, TechnicalReport CS-93-103, Colorado State University (1993); Busetti, Geneticalgorithms overview, http://www.geocities.com/francorbusetti/gaweb.pdf;Beasley et al., An overview of genetic algorithms: Part 1, Fundamentals,University Computing, 15(2), 58-69 (1993); Beasley et al., An Overviewof Genetic Algorithms: Part 2, Research topics, University Computing,University of Cardiff, 15(4), 170-181 (1993); Ingber et al., GeneticAlgorithms and Very Fast Simulated Annealing: Comparison, Mathematicaland Computer Modelling 16(11), 87-100 (1992); Goldberg, GeneticAlgorithms in Search, Optimization, and Machine Learning, Addison-WesleyPub. Co. (1989); Houck et al., A Genetic Algorithm for FunctionOptimization: a Matlab Implementation, NCSU-IE TR 95-09 (1995).

The following is an implementation of the genetic algorithm as itapplies to the present invention. As an example, an optimization problemmay seek the maximal value of the objective function

Val=ƒ(α₁,α₂,α₃)  (27)

This objective function itself is an appropriate fitness functionbecause its value measures how good the solution is. The parameters inthe objective function are consequently encoded to the chromosome asgenes. It is believed a real number representation is more efficientcompared to the binary version. The chromosome with a number ofreal-numbered genes is given as

{α₁,α₂,α₃}  (28)

The first step in genetic algorithm is to give an initial populationwith a number of individuals (that is, with different chromosomes) toreflect the variety. This population is usually generated randomly. Thisstep in fact means to randomly sample the parameter space. Differentfrom random searching technique, the genetic algorithm realizes implicitparallel computation so that it has a faster convergence rate. In termsof schema theory, the fitness of an individual not only represents itsoptimality at that specific position but gives fitness of someunderlying pattern known as schema. After the population has beeninitialized, the fitness score of each individual is evaluated tofacilitate the reproduction step.

Reproduction reflects the natural reality that better fitted individualshave higher probabilities to appear in the next generation. Given thefitness scores of individuals, the next generation is selected in termsof several selection rules. A common method is to select the nextgeneration through the “Roulette Wheel” rule. If the individual i hasthe fitness score Val, then its probability to be selected into nextgeneration is given as

$\begin{matrix}{P = \frac{{Val}_{i}}{\overset{\_}{V}{al}}} & (29)\end{matrix}$

where Val is the average fitness of the population with M individuals

$\begin{matrix}{{\overset{\_}{V}{al}} = \frac{\sum\limits_{i = 1}^{M}\; {Val}_{i}}{M}} & (30)\end{matrix}$

The selection process can be viewed as putting the individuals on awheel and each individual occupies a pie representing its probability.Randomly selecting M individuals from the wheel results in the nextgeneration.

Crossover is the process that one offspring is generated by two parents.The offspring could become a worse-fitted individual. However, there isa good chance that the offspring has better fitness than both of itsparents. Crossover is believed to be the major operator that makes thegenetic algorithm have the capability to investigate the parameter spaceefficiently. One special and efficient real crossover operator is themathematical crossover and is defined as

{α₁ ¹,α₂ ¹,α₃ ¹ }

{p·α ₁ ¹+(1−p)·α₁ ² ,p·α₂ ¹+(1−p)·α₂ ², p·α₃ ¹+(1−p).α₃ ²}{α₁ ²,α₂ ²,α₃²}

{(1−p)·α₁ ¹ ,+p·α₁ ², (1−p).α₂ ¹ +p·α₂ ², (1−p).α₃ ¹ +p·α₃ ²}  (31)

where p is a uniform random number from zero to one.

Mutation is another operator to change the chromosome of an individual.It gives each gene in the chromosome a small probability to be set to arandom number. This operator gives any region in the parameter spacesome probability to be investigated. It therefore helps to escape thelocal minima. However, over-emphasis on the mutation operator leads to astochastic random search.

Starting from the initial population and going over one generation afteranother, finally the genetic algorithm converges to the optimalsolution. As to the optimization problem occurred in shape-basedreconstruction, the objective function is selected as the fitnessfunction. Its parameters is coded into chromosome as

$\begin{matrix}\begin{Bmatrix}{{f_{d\; 1}\left( {X_{1},Y_{1}} \right)} - {f_{d\; 1}\left( {X_{2},Y_{2}} \right)} - \ldots - {f_{d\; 1}\left( {X_{N \times N},Y_{N \times N}} \right)} - {{mp}\; 1} -} \\{{z\left( {X_{1},Y_{1}} \right)} - {z\left( {X_{2},Y_{2}} \right)} - \ldots - {z\left( {X_{N \times N},Y_{N \times N}} \right)} - {{bp}\; 2}}\end{Bmatrix} & (32)\end{matrix}$

Reasonable constraints are added to the parameters in preparation of theoptimization algorithm. First, the region of support of the lesion inx-y plane is readily available in terms of the epi-illuminationnevoscope image (FIG. 4A). This further reduces the number of parametersto represent a surface from N x N to a smaller set. As a consequence,the optimization algorithm has a faster convergence rate. Second, theblood region typically comprises thin layers within a few hundreds ofmicrometers which put a constraint on z(X, Y) . Third, the fractions ofmelanin and blood are not free parameters. They can also be boundedaccording to the appearance of melanoma and the clinical experience.Last, multi-spectral imaging provides implicit constraints. Due to thedistinct absorption spectra of blood and melanin, a reasonable solutionmust also satisfy the measurements of all involved wavelengths.

In summary, the method of solving the shape-based reconstruction problemusing genetic algorithm to involves the following steps:

Step 1: Initialize the population. The values of genes are randomlygenerated as well as subject to some specific constraints. The size ofthe population is sufficient large to reflect randomness.

Step 2: Evaluate the fitness. Each chromosome in the population isevaluated in terms of the objective function and, as a result,associated with a fitness score. This step implicitly includes samplingthe shape-representation to fill into the forward model.

Step 3: Evolve the population. This step includes reproduction in termsof the fitness score, arithmetic crossover and mutation.

Step 4: Check the termination condition. The termination condition inthe experiment is the number of iterations of the algorithm. Once itexceeds the pre-defined number, the algorithm stops. Otherwise, itreturns to step 2 and continues to evolve.

Experiments/Results

To validate the shape-based multi-spectral algorithm, a double-surfacemodel was created to represent malignant melanoma. The first and secondsurfaces are described by a mixed Gaussian function which are given as

ƒ₁(x,y)=MAX(peak1□G(x,y,μ_(1α),μ_(2α),σ_(α)),peak2□G(x,y,μ_(1b),μ_(2b),σ_(b)))

ƒ₂(x,y)=MAX (peak3□G(x,y,μ_(1α),μ_(2α),σ_(α)),peak4□G(x,y,μ_(1b),μ_(2b),σ_(b)))  (33)

where the Gaussian function is

$\begin{matrix}{G = {\left( {x,y,\mu_{1},\mu_{2},\sigma} \right) = {\frac{1}{2{\pi\sigma}^{2}}{\exp \left( \frac{\left( {x - \mu_{1}} \right)^{2} + \left( {y - \mu_{2}} \right)^{2}}{2\sigma^{2}} \right)}}}} & (34)\end{matrix}$

The parameters used in Equation (33) are

$\begin{matrix}\left\{ \begin{matrix}{\mu_{1a} = {0.0375 \times 2\mspace{14mu} {cm}}} & {\mspace{14mu} {{- \mu_{2a}} = {0.0375 \times 3\mspace{14mu} {cm}}}} \\{{- \mu_{1b}} = {0.0375 \times 2\mspace{14mu} {cm}}} & {\mspace{14mu} {\mu_{2b} = {0.0375 \times 3\mspace{14mu} {cm}}}} \\{{\sigma_{a} = {0.0375 \times 2.5\mspace{14mu} {cm}}}\mspace{20mu}} & {\sigma_{b} = {0.0375 \times 2.5\mspace{14mu} {cm}}} \\{{{peak}\; 1} = {100 \times 6\mspace{14mu} {µm}}} & {\mspace{14mu} {{{peak}\; 2} = {100 \times 4\mspace{20mu} {µm}}}} \\{{{peak}\; 3} = {100 \times 8\mspace{14mu} {µm}}} & {\mspace{14mu} {{{peak}\; 4} = {100 \times 6\mspace{20mu} {µm}}}}\end{matrix} \right. & (35)\end{matrix}$

The fraction of melanin is set to 5% between the epidermal layer and thefirst surface ƒ₁(x,y). The fraction of blood is set to 20% between thefirst surface ƒ₁(x,y) and the second surface ƒ₂(x,y). This model hassufficient variation in order to verify the reliability of thereconstruction algorithm. FIG. 8A displays the 3D view of this model.

The double-surface continuous model is sampled and two measurementsΔM⁵⁸⁰ and ΔM⁸⁰⁰ are calculated by Monte Carlo simulation at 580 nm and800 nm respectively. Next a 9×9 rectangular grid is overlapped on theepidermal layer. The region of support of the lesion is counted as 10discrete control points. As a result, the chromosome contains 22 genesand is coded as

$\begin{matrix}\begin{Bmatrix}{{f_{d\; 1}\left( {X_{1},Y_{1}} \right)} - {f_{d\; 1}\left( {X_{2},Y_{2}} \right)} - \ldots - {f_{d\; 1}\left( {X_{10},Y_{10}} \right)} - {{mp}\; 1} -} \\{{z\left( {X_{1},Y_{1}} \right)} - {z\left( {X_{2},Y_{2}} \right)} - \ldots - {z\left( {X_{10},Y_{10}} \right)} - {{bp}\; 2}}\end{Bmatrix} & (36)\end{matrix}$

The fitness function of the genetic algorithm is given as

F _(obj)=α₁ □F _(obj) ⁵⁸⁰ α₂ □F _(obj) ⁸⁰⁰⁺  (37)

where α₁ and α₂ can be chosen on the relative weight assigned to melanincontents to the entire volume. About 30% of the entire volume to bepigmented and used α₁=0.3 and α₂=1 was considered in the experiment.

Four simulations with different constraints were implemented in realnumber presented genetic algorithm. Results are summarized in Table 1.The recovered surfaces are displayed in FIG. 8B-E. In each case, theleft is the first surface and the right is the second surface. There isno constraint for the first surface while the thickness between thefirst surface and the second surface is set to be 300 μm to represent athin layer of blood net.

In an example there are 100 individuals in a population and the geneticalgorithm terminates after 150 iterations. Now referring to FIG. 9, thefitness scores are normalized and displayed with respect to the numberof iteration. As can be seen the genetic algorithm converges in allcases. Now referring to FIGS. 10A and 10B, the deformation process isshown, suggesting that the genetic algorithm has successfully recoveredboth surfaces.

To further evaluate the reconstruction result, volume deviations Volerr1and Volerr2 were introduced. They are defined as

$\begin{matrix}{{{Volerr}\; 1} = \frac{{{{vol}\; 1_{c}} - {{vol}\; 1_{m}}}}{{vol}\; 1_{m}}} & (38)\end{matrix}$

where vol1 _(c) is the calculated volume bounded by the first surfaceand vol1 _(m) is the corresponding volume from the model and

$\begin{matrix}{{{Volerr}\; 2} = \frac{{{{vol}\; {2\;}_{c}} - {{vol}\; 2_{m}}}}{{vol}\; 2_{m}}} & (39)\end{matrix}$

where, vol2 _(c) is the calculated volume bounded by the first andsecond surfaces and vol2 _(m) is the corresponding volume from themodel.

TABLE 1 Summary of Reconstruction Results Melanin Blood RecoveredRecovered Bounds Bounds Melanin blood Volerr1 Volerr2 Case (%) (%) (%)(%) (%) (%) (b) 5-5 20-20 5 20 2.68 16.58 (c) 4.5-5.5 10-30 5.10 18.333.81 20.89 (d) 4-6 10-30 4.64 18.97 5.18 24.71 (e) 3-7 10-30 3.37 10.0844.09 63.50

As shown in Table 1, the constraints have significant impacts on thereconstructed surfaces. Except the loosest constrained case (e), allcases present reasonable reconstructions which are consistent with themodel. Moreover, the reconstructed first surface has a smaller volumeerror than the second surface. There are several reasons to explain thelarger volume error of the second surface. First, the absorptioncoefficient of blood is smaller than that of melanin. As a result, thechange in blood region has less contribution to the fitness function.Second, since a reflectance geometry is adopted in the nevoscope, thesensitivity decreases at deeper layers. This also influences theaccurate reconstruction of the second surface. Third, because the twosurfaces are attached together, the error resulting from the firstsurface would inevitably propagate to the second surface.

The inventive reconstruction techniques are shown to be a good solutionto the ill-posed inverse problem. Using additional wavelengths increasesthe ability to isolate the true solution from a number of possiblesolutions.

Monte-Carlo Simulation-Based Mapping of Multispectral Image Informationonto a 3-D Volume.

Now referring to FIG. 11, in another embodiment a Monte-Carlosimulation-based mapping of multispectral image information onto a 3-Dvolume is provided. The multispectral image information obtained can beprojected onto a 3-D discrete reconstruction volume through theknowledge-based model derived from a Monte-Carlo simulation providingprobability distribution functions of energy and path length and energyin the 3-D sub-surface space. The Monte-Carlo simulation of theskin-tissue model (without the lesion) provides a probabilitydistribution of each photon path between the source and detector in thesimulated medium. An iterative ART algorithm can be used to employ thismatrix of initial weights and path-length information to compute thediffused reflectance distribution that is compared with the acquiredmultispectral images of skin of the human subjects (near theskin-lesion). Thus, the weights and path-length information of thereconstruction space of the skin is updated using an AlgebraicReconstruction Technique (ART) (See, S. Srinath Maganti,“Three-Dimensional Optical Tomographic Image Reconstruction fromNevoscope Images”, Ph.D. Dissertation, Advisor: Atam P. Dhawan,University of Cincinnati (1997); S. Maganti and A. P. Dhawan, “3-DNevoscope image reconstruction using diverging ray ART”, ProceedingsSPIE International Conference on Biomedical Optics (1997)) forminimization of the error between the computed diffused reflectance andacquired multispectral images of skin. This updated weight matrix andpath-length information in the 3-D reconstruction space now serves asinitial weight matrix for reconstructing the skin-lesion informationfrom multispectral images of the skin-lesion using the iterative ART.Thus, Monte Carlo simulation provides the forward model of imaging forthe voxel weights and path-length information which is then used inreconstruction through an iterative Algebraic Reconstruction Technique(ART) (See, S. Srinath Maganti, “Three-Dimensional Optical TomographicImage Reconstruction from Nevoscope Images”, Ph.D. Dissertation,Advisor: Atam P. Dhawan, University of Cincinnati (1997); S. Maganti andA. P. Dhawan, “3-D Nevoscope image reconstruction using diverging rayART”, Proceedings SPIE International Conference on Biomedical Optics(1997)) as depicted in FIG. 11.

EXAMPLE

It is assumed a skin lesion property (such as melanin, hemoglobin, orde-oxy-hemoglobin concentration) is represented by F_(i)(x,y,z); i=1, .. . , I where I is the number of skin-lesion properties. It is alsoassumed that the Monte-Carlo simulation and iterative ART-basedreconstruction of skin provides the updated initial weight matrix with3-D distributions of photon energy density represented by E_(j)(α,β,γ)with path lengths P_(k,i,j)((α₁β_(m),γ_(n))⁰, E_(j;1mn)),(α₁,β_(m),γ_(n))¹, E_(j;lmn)), . . . , (α₁,β_(m),γ_(n))^(Z),E_(j;1mn))); k=1, . . . , K where K is the number of paths betweensource and detector for a specific imaging wavelength j (j=1, . . . , J)for imaging the skin-lesion property F_(i)(x,y,z). Each path length is alist of voxel coordinates (α₁,β_(m),y_(n))^(Z); z=1, . . . , Z (where Zis the number of voxels between source and detector) with theirrespective associated energy level E_(j;lmn). The normalized probabilitydistribution of photon energy provides initial weights forreconstructing the object information F_(i)(x,y,z) as RF_(i)(x′,y′,z′)where RF_(i) is the respective reconstructed property in the 3-Dreconstruction space (x′,y′,z′). Let NEVO₁₃ TRANS function describe theimaging the skin-lesion property, F_(i)(x,y,z) through NevoscopeTransillumination with the forward model:

[M_(j)(α,β,γ)]=NEVO_TRANS{F_(i)(x,y,z)}=[F_(i)(x,y,z)].[W_(j)(α, β,γ)]where [M, (α,β,γ)] is the measured or acquired image, [W_(j)(α,β,γ)] isthe updated weight matrix, which is approximated and modeled throughMonte-Carlo simulation with [E_(j)(α,β,γ)].

The objective of the reconstruction process is to obtain an estimateRF_(i)(x′,y′,z′) of the object property F_(i)(x,y,z). The ART algorithmstarts projecting the acquired image data onto the reconstruction spaceof voxels using the initial weights and path length information obtainedfrom the Monte-Carlo simulation (forward model). In each iteration, thevoxel weight and path length information is used to assign new weightvalues such that the error between the reconstructed diffusedreflectance output from the medium using the new voxel weights and theacquired image is minimized. Thus, the reconstruction process concludesthat the incident photon energy radiation on the reconstructed objectprovides the diffused reflectance output as close as possible to theacquired image as

[RF _(j)(x′,y′,z′)]=[P _(k,i,j)((α₁,β_(m),γ_(n))⁰ ,E _(j;1mn)), (α₁,β_(m),γ_(n))¹,E _(j;l,mn)), . . . , (α₁,β_(m),γ_(n))^(Z) ,E_(j;lmn)))]^(INV).

[M _(j)(α,β,65 )]=[W _(j)(α,β,γ)]^(INV).[M _(j)(α,β, γ)]

Thus, through the ART reconstruction process an inverse of the weightmatrix, [W_(j)(α,β,γ)]^(INV) is obtained through the iterative process.However, since the matrix is very large to invert, the convergence canbe obtained using a hierarchical partitioning. Also, the a prioriknowledge of the path length from the Monte-Carlo simulation can be usedto reduce the size of the sub-matrices to be inverted. Thus, theassignment of the acquired image data onto a reconstructed 3-D space canbe efficiently performed for each property of the skin-lesion throughspecific wavelength based imaging. Thus, the individual wavelength basedreconstruction from each of the acquired multispectral images will beobtained. The specific property—(such as melanin, hemoglobin, etc.)based 3-D reconstructed distributions would be obtained throughcombining the multispectral reconstructions weighted by their respectivevolume fraction model described in the equations above. The totalmelanin fraction would be obtained from multispectral images as

$\begin{matrix}{\left\lbrack {{RF}_{melanin}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)} \right\rbrack = {{\sum\limits_{j = 1}^{J}\; {{\left( {f_{j,{epmel}} + f_{j,{pmel}} + f_{j,{rmel}}} \right)\left\lbrack {{RF}_{j}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)} \right\rbrack}\left\lbrack {{RF}_{{oxy}\text{-}{hemo}}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)} \right\rbrack}} = {{\sum\limits_{j = 1}^{J}\; {{\left( {f_{j,{poxy}} + f_{j,{roxy}}} \right)\left\lbrack {{RF}_{j}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)} \right\rbrack}\left\lbrack {{RF}_{{de}\text{-}{oxy}\text{-}{hemo}}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)} \right\rbrack}} = {\sum\limits_{j = 1}^{J}\; {\left( {f_{j,{pdeoxy}} + f_{j,{rdeoxy}}} \right)\left\lbrack {{RF}_{j}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)} \right\rbrack}}}}} & \;\end{matrix}$

where f-_(j,epmel, f) _(j,pmel), and f_(j,rmel), are the wavelengthdependent volume fractions of melanin in, respectively, epidermis,papillary dermis and reticular dermis; and f_(j,poxy), f_(j,roxy),f_(j,pdeoxy), and f_(j,rdeoxy) are the wavelength dependent volumefractions of oxy-hemoglobin and de-oxyhemoglobin in, respectively,papillary dermis and reticular dermis. Validation of the melanin andblood-flow volume information can be compared with the pathology fortest cases for validation.

Experiments

Different visible light and near-infrared wavelengths in the range of400 nm 1000 nm (preferably 500 nm to 960 nm) may be used to obtaindiffused reflectance-based absorption measurements with different levelsof depth of penetration to estimate volumes of selected chromophoressuch as melanin, oxy-hemoglobin and de-oxy-hemoglobin. A summary of thedepth of penetration and characteristic absorption to chromorphores ofinterest is provided in Table 2.

TABLE 2 Specific light wavelengths and their properties for use inmultispectral skin imaging. Wavelength Optical Properties in Skin 405 nmSuperficial penetration (<0.5 mm), used in fluorescence excitation ofALA. High melanin absorption. 470 nm Superficial penetration (<1.0 mm),used in fluorescence excitation of fluorescein, highly absorbed inmelanin. Excellent for imaging superficial melanin associated with sundamage. 500 nm Superficial penetration (<1.2 mm), used for fluorescenceimaging 520 nm Some fluorescence imaging (<1.5 mm penetration), 560 nmSuperficial penetration (<1.8 mm), highly absorbed by de-oxy-hemoglobin.Highlights superficial blood vessels. 580 nm Deeper penetration intissue (<2.3 mm), excellent for imaging most pigmented lesions withoutinterference from the superficial melanin, peak for oxy- hemoglobinabsorption. 610 nm Deeper penetration in tissue (<2.5 mm), excellent forcontrast imaging deoxygenated hemoglobin and veins due to largedifference in absorption for oxy and de-oxy-hemoglobin. 660 nm Deeppenetration in tissue (>3.0 mm), shows deep pigmentation; largedifference in absorption for oxy and de-oxy-hemoglobin. 760 nm Very deeppenetration in tissue, peak for de-oxy-hemoglobin and water absorption.800 nm Very Deep penetration, almost equal absorption for de-oxy andoxy-hemoglobin. 910 nm Very Deep penetration in tissue, peak foroxy-hemoglobin absorption.

Optical wavelength of 560 nm provides superficial penetration (<1.8 mm)with high absorption by deoxy-hemoglobin and is suitable for imagingsuperficial blood vessels. While 610 nm wavelength provide deeperpenetration in tissue (<2.5 mm), it is also excellent for contrastimaging deoxygenated hemoglobin and veins due to large difference inabsorption for oxy and de-oxy-hemoglobin. The images obtained from 660nm and 910 nm provide measurements for the blood oxygen saturation level(SO₂). Thus a combination of selected wavelengths from spectralmeasurements of specific chromophores can provide useful informationabout melanin, oxy- and deoxy-hemoglobin.

For measurements related to melanin, vascular absorption due tode-oxyhemoglobin, and blood oxygen saturation level, 560 nm, 580 nm, 610nm, 660 nm, 760 nm, 800 nm and 910 nm wavelengths are selected forimaging through a multispectral nevoscope apparatus 200. The imagesobtained from 660 nm, 760 nm, 800 and 910 nm can be used to estimate theblood oxygen saturation level (SO₂) of the skin-lesion. Thesemeasurements can be used in 3-D reconstruction of Hb and HbO₂distributions using a Monte-Carlo simulation based reconstruction methodas described herein. Other wavelengths (560 nm, 580 nm, and 610 nm) canbe used for visual and 2D analysis of melanin and superficial-to-deepdeoxy-hemoglobin-based vascular absorption.

Near-Infrared Imaging for Oxygen Saturation Level Measurement (HbO₂ andHb absorption)

In the absorption spectra of HbO₂ and Hb, as the oxygen saturation goesup, more of the light longer than 800 nm gets absorbed while the shorterwavelength light absorption decreases. To compare the localizedabsorption difference, we need to obtain the light spectrum with aspecific optical path length, and accurately measure how much spectraldistortion is caused by the blood absorption, for long and shortwavelength range respectfully. In order to analyze the wavelengthdependent absorption coefficients, the image obtained from the CCDcamera for a fixed wavelength will be analyzed to estimate theabsorption coefficients at various wavelength.

Generally the light signal level, r_(s) decays with increasing depth.The attenuation coefficient is the summation of scattering coefficientα_(s) and absorption coefficient α_(α)(k). To simplify this problem, weassume α_(s) is constant over the bandwidth of the light used. Thecumulative spectral distortion at depth l is determined by the lightabsorption from the end of probe to l, which can be written as:

{circumflex over (r)}_(s)(k,l)=r _(s0)exp {−∫₀ ^(l)[α_(s)(l′)+α_(α()k,l′)]dl′}

By averaging {circumflex over (r)}_(s)(k,l) in long and short wavelengthrange, the following are obtained:

${\hat{r}}_{s\_ {short}}\left( {{l} = {r_{s\; 0}{\exp\left( {2{\int_{0}^{l}{{\alpha_{s}\left( {l^{\prime}{l^{\prime}}} \right)}\left\{ \overset{\_}{\exp\left\lbrack {\int_{0}^{l}{\alpha_{a}\left( {k{l^{\prime}}} \right\rbrack}} \right.} \right\}_{short}\mspace{20mu} r_{s\; 0}{\exp \left( {\int_{0}^{l}{{\alpha_{s}\left( l^{\prime} \right)}{l^{\prime}}}} \right)}{A_{short}(l)}{{\hat{r}}_{s\_ {long}}\left( {{l} = {r_{s\; 0}{\exp \left( {2{\int_{0}^{l}{{\alpha_{s}\left( l^{\prime} \right)}l^{\prime}}}} \right)}\left\{ \overset{\_}{\exp\left\lbrack {\int_{0}^{l}{\alpha_{a}\left( {k{l^{\prime}}} \right\rbrack}} \right.} \right\}_{long}\mspace{20mu} r_{s\; 0}{\exp \left( {\int_{0}^{l}{{\alpha_{s}\left( l^{\prime} \right)}{l^{\prime}}}} \right)}{A_{long}(l)}}} \right.}}}} \right.}}} \right.$

{⁻}_(short,long) indicates that the average value of the quantity inside of the brackets in short (e.g. 660 nm) or long wavelength range(e.g. 910 nm).

Although the spatial dependence of absorption coefficient is not exactlyknown, A_(short) and A_(short) are related to the oxygen saturationlevel, because the blood absorption coefficient α_(a)(k,l) is a functionof SO₂

α_(a)(k,l)=[SO ₂(l)]α_(HbO) ₂ (k)+[1−SO ₂(l)]α_(Hb)(k)

where α_(HbO2) is the wavelength dependent absorption coefficient ofHb)₂ corresponding to wavevector k; α_(Hb) is the absorption coefficientof Hb. And SO₂ varies as l. R(l) defined as difference of absorptionbetween long and short wavelength light can be computed as

${R(l)} = {\frac{{{\hat{r}}_{s\_ {short}}(l)} - {{\hat{r}}_{s\_ {long}}(l)}}{{{\hat{r}}_{s\_ {short}}(l)} + {{\hat{r}}_{s\_ {long}}(l)}} = \frac{{A_{short}(l)} - {A_{l{ong}}(l)}}{{A_{l{ong}}(l)} + {A_{short}(l)}}}$

R(l) is not determined by the absolute amplitude of {circumflex over(r)}_(s)(k,l), but determined by the difference in absorption of thelong and the short wavelength range. Therefore, R(l) does not show howmuch spectral absorption is caused by blood, but rather shows the SO₂related difference in spectral absorption. By images obtained from theshort and long wavelength, processing data with aforementionedalgorithm, we can obtain 2D SO₂ image with R(l_(x),l_(z)).

Computer Classification of Multispectral Nevoscope Images

The disclosed classification system uses an adaptive fuzzy c-meansclustering (AFCM) algorithm to provide an optimal set of classifiableclusters of the correlated features extracted from the white-light andmultispectral nevoscope images. Now referring to FIG. 12 a schematicdiagram of a computer classification system using adaptive Fuzzyclustering and portioning is detailed.

Adaptive Fuzzy c-means Clustering Method

The fuzzy c-means clustering method utilizes an adaptable membershipvalue that can be updated based using the distribution statistics of thedata points assigned to the cluster minimizing the following objectivefunction J_(m)(U,v)

${J_{m}\left( {U,v} \right)} = {{\sum\limits_{i = 1}^{c}\; {\sum\limits_{j = 1}^{n}\; {u_{ij}^{m}d_{ij}^{2}}}} = {\sum\limits_{i = 1}^{c}\; {\sum\limits_{j = 1}^{n}{u_{ij}^{m}{{x_{j} - v_{i}}}^{2}}}}}$

where c is the number of clusters, n the number of data feature vectors,u_(ij) is the fuzzy membership, m is the fuzziness index, x_(j); j=1 . .. , n is the input data feature vectors, v_(i); i=1, . . . c are thecluster centroids; and d_(ij) is the Euclidean distance vector. Based onthe constraints defined on the distribution statistics of the datapoints in the clusters, fuzziness index can be defined between 1 and avery large value (maximum allowable variance within a cluster). Themembership values can be defined as:

0 ≤ u_(ij) ≤ 1  for  all  i, j${\sum\limits_{i = 1}^{c}{u_{ij}\mspace{14mu} {for}\mspace{14mu} {all}\mspace{14mu} j\mspace{14mu} {and}\mspace{14mu} 0}} < {\sum\limits_{j = 1}^{n}u_{ij}} < {n\mspace{14mu} {for}\mspace{14mu} {all}\mspace{14mu} i}$

An adaptive version of the FCM algorithm using a weighted derivativeinformation is described by Pham and Prince, D. L. Pham, J. L. Prince,“Adaptive fuzzy segmentation of magnetic resonance images,” IEEE Transon Med. Imaging, 18, 9, 737-752, (1999). The adaptive FCM algorithmprovides a good regularization in the clustering process through thefirst- and second-order derivative terms that force the bias field to besmooth. As the objective function is minimized for clustering, theminimization of the derivative information provides smooth homogenousclusters.

The objective function J_(AFCM)(U, v) to be minimized for an adaptivefuzzy c-means clustering algorithm can be given as

${J_{ACFM}\left( {U,v} \right)} = {{\sum\limits_{i = 1}^{c}{\sum\limits_{j = 1}^{n}{u_{ij}^{m}{{x_{j} - {g_{j}v_{i}}}}^{2}}}} + {\lambda_{1}{\sum\limits_{j = 1}^{n}{\sum\limits_{r = 1}^{R}\left( {D_{r} \otimes g} \right)_{j}^{2}}}} + {\lambda_{2}{\sum\limits_{j = 1}^{n}{\sum\limits_{r = 1}^{R}{\sum\limits_{s = 1}^{R}\left( {D_{r} \otimes D_{s\;} \otimes g} \right)_{j}^{2}}}}}}$

where g_(j) is the unknown gain field to be estimated during theiterative clustering process, m is the fuzzification index (m>1), R isthe dimension of the image, D_(r) is the derivative operator along ther^(th) dimension of the image, and λ₁ and λ₂ are, respectively, theweights for the first-order and second-order derivative terms. Secondand third terms in the objective function provide the smoothnessconstraints.

It should be noted that though there may be as many as 14 classes forthe lesion classification, the number of clusters c may be much largerthan 14 depending on the distribution of data in the feature space. Allclusters are required to be homogenous with respect to each category andmust satisfy the maximum allowable variance criterion. If a clusterviolates the criterion when mapped with the pathology classificationlabels, it is split and a new value of c is computed. After the mappingand labeling of each cluster is completed, the class models aredeveloped based on the set of clusters and their parameters includingmean, variance and fuzzy membership functions. For example, atwo-dimensional feature space is shown in FIG. 13A for simplicity toillustrate a hypothetical class distribution. Class A includes twoclusters, 1 and 3, while the Class C includes clusters 4 and 6. Themodel for respective classes must include all clusters and theirparameters. These models are stored in the database for each class andupdated on a need basis whenever there is disagreement between thepathology and the computed classification from any case in the test set.Thus, the model adapts to the new knowledge and learns new modelparameters re-organizing the clusters. The new learning is achievedthrough the competitive learning approach using a leader-followerclustering algorithm, see, e.g., A. K. Jain and R. C. Dubes, Algorithmsfor Clustering Data, Prentice Hall, Englewood Cliffs, N.J., (1998). Thenearest clusters are polled for the fuzzy membership functions. Based onthe “Winner-Take-All” strategy, the cluster with maximum fuzzymembership function is chosen for the inclusion of the new data vectoras long as it meets the cluster homogeneity and variance criteria.

After the models for pathology classification classes are constructedthrough labeling of corresponding clusters and their parameters, anintegrated competitive learning approach for the classification isdeveloped. In one example hyperplanes as shown can be placed by straightlines in FIG. 13A. Hyperplanes can be placed in the feature spacethrough checking the clusters in each model and creating convex setsaround them. W. Grohman and A. P. Dhawan, “Fuzzy convex set basedpattern classification of mammographic microcalcifications”, PatternRecognition, vol. 34(7), pp. 119-132 (2001). The creation of convex setsis time consuming but if the cluster parameters are stored for eachmodel class, their respective fuzzy membership functions can be used topartition the feature space efficiently. After the convex sets aroundthe clusters are identified, hyperplanes can be placed for computing thefuzzy classification membership functions to allow final classificationusing a “Winner-Take-All” strategy as shown in FIG. 13B.

Classification

A fuzzy classification membership function M_(f) would be constructed toreflect the true shape of the convex subset with the greatest precisionpossible. The function M_(ƒ) is defined for each subset ƒ(ƒf=1, 2, . . ., K) of the class model as follows:

${M_{f}(x)} = {{\sqrt[L_{f}]{{\prod\limits_{i = 1}^{L_{f}}\; \theta_{i}},}\mspace{14mu} \theta_{i}} = \frac{1}{\left( {1 + ^{\lambda_{if}\phi_{i}x}} \right)}}$

where: L_(ƒ) is the number of separating hyperplanes for the subset ƒ,φ_(i) is the i-th separating hyperplane function for the subset, in thevector form, x is the input vector in the augmented form, and λ_(iƒ) isthe steepness (scaling) coefficient for the i-th hyperplane in thesubset ƒ. The value of λ_(iƒ) depends on the depth of convex subset ƒ,as projected onto the separating hyperplane H_(i) (defined by φ_(i)):

${\lambda_{if} = \frac{- {\log \left( \frac{1 - \chi}{\chi} \right)}}{\mu_{if}}},{u_{if} = {\frac{1}{N}{\sum\limits_{j = 1}^{n}{\phi_{1}x_{j}}}}}$

where: n is the number of training points in the covex subset ƒ, φ_(i)is the i-th hyperplane equation in the vector form, μ_(iƒ) is the depthof the convex subset f, as projected onto i-th hyperplane, x_(j) is anaugmented coordinate vector of the j-th point in the subset, and x isthe center value of the membership function.

The output CO of the classifier is the category C of the convex setfuzzy classification membership function M, that attains the highestvalue for the specified input pattern x, i.e.:

$\left( {{CO} = {C\begin{matrix}{1 \leq \overset{\forall}{f} \leq K_{{{M_{f}{(x)}} < {M_{i}{(x)}}},{M_{i} \in C}}} \\{f \neq i}\end{matrix}}} \right)$

where: CO is the output of the classifier, K is the number of convexsets obtained during training (number of fuzzy function neurons in thefuzzy membership function layer), M_(i) is the highest fuzzyclassification membership function value for the input x, and C is themodel class of convex subset used to construct membership functionM_(i).

Applicants have attempted to disclose all embodiments and applicationsof the described subject matter that could be reasonably foreseen.However, there may be unforeseeable, insubstantial modifications thatremain as equivalents. While the present invention has been described inconjunction with specific, exemplary embodiments thereof, it is evidentthat many alterations, modifications, and variations will be apparent tothose skilled in the art in light of the foregoing description withoutdeparting from the spirit or scope of the present disclosure.Accordingly, the present disclosure is intended to embrace all suchalterations, modifications, and variations of the above detaileddescription.

All references cited herein are incorporated fully by reference.

What is claimed is:
 1. A multispectral nevoscope for examining in situ askin lesion of a patient, comprising: (a) a housing having a centralaxis and a base with a central aperture for positioning over a portionof skin of a patient, the portion of skin including the lesion to beexamined; (b) a face plate positioned in register with the centralaperture comprising an imaging area having plural illumination fiberchannels operably connected to an illuminator-filter assembly and pluralreceiver fiber channels; and (c) at least one transillumination devicedisposed proximal to the faceplate operable to transmit light in aconverging ring into the skin of the patient in the region of the lesionto uniformly transilluminate the region of skin of the patient includingthe lesion.
 2. The multispectral nevoscope of claim 1 wherein thetransillumination device comprises at least one ring illuminatordisposed on the face plate.
 3. The multispectral nevoscope of claim 2wherein the at least one ring illuminator comprises a first outer ringilluminator comprising at least one transillumination fiber channelhaving a first wavelength oriented at 45 degrees to produce a convergentbeam for multispectral or specific wavelength-based background imageexcitation, and a second outer ring illuminator comprising at least onetransillumination fiber channel having a second wavelength oriented at45 degrees to produce a convergent beam for multispectral or specificwavelength-based background image excitation.
 4. The multispectralnevoscope of claim 3 wherein the first wavelength is shorter than thesecond wavelength.
 5. The multispectral nevoscope of claim 2 furthercomprising a source mask wherein the imaging area illumination fiberchannels and receiver fiber channels are in positions in a matrixcorresponding to the source mask.
 6. The multispectral nevoscope ofclaim 1 further comprising a focusing lens and optionally a polarizer.7. The multispectral nevoscope of claim 1 further comprising amultispectral filter.
 8. The multispectral nevoscope of claim 1 whereinthe face plate is removably attached to the housing and is operablyconnected to a removable lens and/or polarizer wherein the lens and/orpolarizer provides an interface to the imaging area.
 9. Themultispectral nevoscope of claim 1 wherein at least one of the pluralreceiver fiber channels is operably connected to a CCD image sensor. 10.The multispectral nevoscope of claim 1 comprising a device for recordingimages of the transilluminated skin lesion.
 11. The multispectralnevoscope of claim 10 comprising a mask to receive the light from fibersnot used for illumination, a focusing lens and/or polarizer operablylinked to the mask and a CCD sensor to record images.
 12. A systemcomprising a multispectral nevoscope according to claim 1 comprising acomputing device operable to classify images collected by themultispectral nevoscope.
 13. The system according to claim 12 whereinthe computing device employs software operable to classify images usingan adaptive fuzzy c-means clustering algorithm to provide an optimal setof classifiable clusters of correlated features extracted fromwhite-light and multispectral nevoscope images.